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ABSTRACT 


The Du Fort-Frankel dif c erence scheme is generalized to difference 
operators of arbitrary high order accuracy in space and to 
arbitrary order of the parabolic differential operator. Spectral 
methods can also be used to approximate the spatial part of the 
differential operator. The scheme is explicit, and it is uncondi- 
tionally stable for the initial value problem. Stable boundary 
conditions are given for two different fourth order accurate 


space approximations. 
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I. INTRODUCTION 


The Du Fort-Frankel difference scheme for solving parabolic 
equations, see e.g. [6], has the advantage of being explicit and 
yet unconditionally stable. The consistency requires that At goes 
to zero faster than Ax does, but this requirement is in practice 
not too severe if the coefficient of the second derivative is 
small. 

An analysis of a semidiscretized parabolic model problem, 
performed along the same lines as in (2), [3] for hyperbolic problems, 
shows that higher order accurate approximations are more efficient 
except for very low requirements on the accuracy of the results. 
Therefore, in this paper the Du Fort-Frankel scheme will be 
generalized to difference operators of arbitrary high order accuracy 
in space and to arbitrary order of the parabolic differential 
operator. The number of space dimensions is also arbitrary and so 
is the number of equations in the system. The scheme is explicit 
and unconditionally stable. For a system with l differential 
equations we also avoid the solution of an £x£ system of equations 
for each gridpoint which would result from a straightforward 
formulation of the scheme. In addition to f inice differences, 
spectral methods and finite element methods can also be used to 
approximate the spatial part of the differential operator in our 
scheme . 

As for the original Du Fort-Frankel scheme, consistency imposes 
a restriction on At in relation to Ax. However, for the type of 
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applications that we have in mind, like the viscous Navier-Stokes 
equations, the dominating truncation error comes from the space dis- 
cretization. Furthermore, when the time dependent equations are used to 
obtain a steady state solution, the truncation error from the time 
discretization is of no importance, assuming that the scheme 
converges for t 

The generalization to higher order accurate approximations in 

space was given by Swarz [8] for the scalar equation u t = ou xx with 

periodic boundary conditions, be studied the efficiency for different 

orders of accuracy and found e.g. that 12th order accurate operators 

-2 

are optimal in a certain sense for a relative precisicr. of 10 , 

and even higher order for higher nrecision. In real applications 
with non periodic boundary conditions, we think that 4th or 6th 
order operators are more realistic. 

In section 2 the scheme will he prese'ted for a r *quence of 
differential equations of increasing generality. In sc-ction 3 the 
stability proofs are given, and in section 4 the so called Fourier 
method is treated. In section 5 the stability of the mixed initial 
boundary value problem is proven for two different 4th order accurate 
space approximations. Section 6 contains a presentation of some 
numerical experiments that were done for the Burger's equation. 


original pam jb 

Of POOR QUALUTJ 
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2. THE DP FORT-FRAHKEL METHOD FOR FINITE DIFFERENCE SCHEMES 

In order to illustrate the idea of the original Du Fort-Frankel 
scheme we start from the simple equation: 


(2.1) u fc ■ ou^ , a > 0. 


Hit*. u" * u(jAx,nAt), it is well known that the scheme 


2 At 


(Ax) 


^7 ( Vr 2u j +u j-i ) 


is unconditionally unstable. However, if we replace u? by 


y(Uj + ^tu” - ^) , the scheme becomes unconditionally stable. For 
higher order approximations to u xx it is not enough to replace u ^ 

by some average. He adopt, therefore, another approach (see also 
Swartz [8]). 

—2 2 a ^ 

Let (Ax) D 2 be a 2p th order approximation to — j ; then the 

3x 


generalized Du Fort-Frankel scheme will be 


( 2 . 2 ) 


2At 


a d 2 „n 


n+l_-j„n. n-1. 


(Ax) 


2 ~2p' 


u 1 ? - -22_ ( u " +1 -2u" + u" 
3 (Ax) * ^ J 3 


where y is to be chosen such that the scheme is unconditionally 

stable. The second term on the right hand side is a stabilizer. 

At 2 

and it is an approximation to y a (^) u fct - 


Therefore, consistency 
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requires that ~ ♦ 0; moreover, in order to minimize the truncation 
error me mould like to choose y as small as possible. It should be 
noted that the operator D* p is not necessarily a difference 
operator; me can use any method of approximation such as spectral 
methods [3] , 151, and finite element methods [7]. For difference 
approximations, y is qiven by 

(2.3) Y > Y 0 = | max |D*(C) | 

Isli* P 

where 


(2.4) D^U) * e~ ikx D* p e ifc: , C - kix . 


In the next section me mill shorn that (2.3) yields unconditional 
stability. The original Du Fort-Frankel ncheme corresponds to 
Y«Yq. However, in this case the stability is not clear for 
variable coefficients. 

We would like to consider in detail some difference approxima- 
tions for D^p, and the first one is the symnetric explicit operator. 


n_ n 


Let D + uV«uy +1 - u y, D - u j“jj* u j-i‘ Then 


P-1 




(2.S) ^2p E D+D- I <-D j «° + D -» 3 

j=0 


From (2.5) we deduce the following formula for the Fourier 


transform 
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*2 P * ' 4 * l " J i 3 | 0 Tf jv jf iijti) Bin23 1 • 


It is obvious that D 2p U) takes its maximum at € 


± tt . Therefore 


(2.6) y, 


P-1 

l 

j*0 


4 j (jl) 2 


(2j+l) 1 (j+1) 


and for every y such that y > y Q , (2.2) is unconditionally stable. 
2 2 

The operators D 2 and are of special interest from the practical 
point of view. Formula (2.5) yields 

(2.7a) D 2 » D + D_ 

(2.7b) D 2 » D + D^(I- D + DJ = Q x 


For the first operator we have Yg = If and for the second one 

4 

Yq * y. Another operator which is of importance is the fourth 
order implicit operator given by the following formula: 


( 2 . 8 ) 


3 2 u 


D D 
+ - 


dx 


(Ax) 


jT 

TT 


n _ 
U D = 




(Ax) 


1 . n 

— ’ °2 u j 


It was shown by Kreiss (5 ) that this operator is more accurate 
than the explicit operator defined in (2.7b). It is easily 
verified that for this case Yq * \ • 



*e consider now a parabolic system of equations of the 


form 


(2.3) » A(x # t,u)u jat , 

where u ,ii u 1 component vector and A is an t*l matrix. 
The condition for uniform parabolic! ty iss 


(2.10) Real X(A) > > 0 

for each eigenvalue X(A) of A. There are two ways to extend the 
method defined by (2.2). The first is: 


u r- u £. 

2At 



— a!?(uI' I>1 -2u?-Ki;“ 1 ) 
(Ax) 2 33 33 


9 


which requires the solution of an system of equations in every 
time step# and therefore, this method is not desired. A better 
method can be obtained by: 


( 2 . 11 ) 


2At 



n 


A j D 2p u j 



p(A®) (Uj* 1 -2u“+u““ 1 ) 


where 

(2.12) p (A) - max |X(A)| , 


and X is given by (2.3). (Here we assume that p(A”) is known 
explicitly, or that a good upper bound is known. Thw spectral 
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radius should not be computed numerically at each timestep.) A 
stability analysis for certain classes of the matrix A will be 
presented in the next section. 

We will now discuss the two dimensional case for systems: 


(2.13) u * A(x,t,u)u + B (x, t,u)u + C(x,t,u)u . 

t AA yy 

The equation is said to be uniformly parabolic in the sense of 
Petrovskii if there exists a constant $ 2 independent of x,t and 
u such that: 


2 2 

(2.14) Real X (Au^+BtOjWj + Cw 2 ) >_ $ 2 > 0 for all real # u> 2 
2 2 

with w^+u) 2 *l. The Du Fort-Frankel scheme for (2.13) is: 


-li ii_ 

2At 


A? fl (D 2 ) u* + 


(Ax) 2 ji 2p x ji (Ay) 


1 _n ,_2 . n 
J c ji ^°2p y u jt 


(2.15) 


B ?o ( 


AxAy “jt'~2p'x'~2p'y ji 


(Ax) 2 


2 ' 
(Ay ) 1 


* (u 


n+1 

ji 


2u jt 
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where (D) x means a difference operator in the x direction and 
(D)y is a difference operator in the y direction. D 2p is any 
approximation to the first derivative accurate up to 2pfch order. 
It should be noted that the stabilizing term, the last term in 

(2.15), is independent of B . This means thac (2.15) is a very 
simple explicit method that can be extended easily to more than 
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two space dimensions. We can determine again v by (2.3). 

The last problem to be treated is a general parabolic 
differential equation of order 2m in s space dimensions: 

(2.16) - Y. A v ( x » t ) ^} 1 ***dg Su » X « (x lt ...,x s ) 

|»i|=2m 

s 

v® lv| = y V t » Sx” 

i=l J 

The equation is said to be parabolic if there is a constant 
such that all the eigenvalues of 

^ AJX.tJft.j)'’ 1 ... (l« s )'’ s 

| vi | -2m 

satisfy 


(2.17) Real X < < 0 

for all « = (<r 1# . . . »® s )» real and Y ®^ = 1 • The scheme 

for (2.16) will be 


(2.18) 
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where In the second sum all the terms with mixed derivatives 
are excluded. 

3. Stability for the Initial value problem . 

In this section it will be shown that the scheme presented 
In the first section is unconditionally stable for the linear 
pure initial value problem. We shall maxe use of the stability 
theory developed by WidlundflOJ and we will assume that the 
reader is familiar with that paper. 

We start with the following lemma: 

Lemma 3.1 . Consider the equation: 

(3.1) X 2 - 1 = -2yX - x(x-l) 2 


where 

(3.2) x > | > 0 . 

Then the roots x + and x_ of (3*1) satisfy 

(3.3) |X ± | < 1 * 

where the equality sign holds only if y - 0, and then only for 

v • 

1 roof : The roots of (3.1) are 

X ± * 1+x ± J 1 " y(2x-y) ] . 
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If l-y(2x-y) > o then by (3.2) we get 

m y(2x-y) < 1 , | x-y | < x , 

and therefore (x^i <,1 • Equality holds in (3*3) only if y=o , 
in this case 



If now 1 - y(2x-y) < o * then and X. are complex and 

therefore 



< 1 . 


This completes the proof. 

We discuss first the method defined by (2.2). The Fourier 
transform of (2.2) is exactly (3*1) with 


y = 




x 


= 2vo 


At 

rr 

(axt 


Equation (3*2) yields the condition (2.3) and by lemma 2.1 only 
one of the roots lies on the unit circle, which is sufficient for 
stability. Moreover, since ^^(rj = 0 only if % » 0 , we 
have also: 


( 3 . 4 ) 


M < i - «l?l 2 • 


The bound on |x| in (3.4) is important if lower order terms 
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are included in (2.1). It should be noted that with v = y q * 
as in the original Du-Fort- Franke 1 scheme, (3.**) is not fulfilled. 

In order to investigate the stability of (2.11) we make 
several assumptions that simplify the analysis. We will assume 
that A is independent of u and has real eigenvalues. The 
second assumption applies to a large class of problems such as 
the compressible viscous Navier StoKes equations. (We recognize, 
though, that the Navier Stos.es equations are not uniformly 
parabolic since the continuity equation does npt contain second 
derivatives) . 

We shall show now that the scheme (2.11) is a parabolic 
difference scheme in the sense of Widlund [10]. Rewrite first, 
the Fourier transform of (2.11) to get: 



Li o j 

The eigenvalues z of G satisfy the equation (3.1) with X = z and 


( 3 . 6 ) 
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and 


(3-7) 


x = o(A) 

(«) 


Again if (2.3) is satisfied we get: 


z) < 1-6|?| , \%\ < ir 


It remains to checK the root condition of Widlund [10], that is 
v;e have to checK that the eigenvalues of 


2x l-x“| 

I+x Tof 


2 0 


are not outside the unit circle and are simple on the unit circle. 

x+i 

But the eigenvalues of H are and therefore, the root 

condition is satisfied. This completes the proof that (2.11) is 
unconditionally stable. 

It should be noted that H defined in (3.8) satisfies a 
stronger condition than the root condition since only one of 
the roots lies on the unit circle. In fact the conditions of 
theorem 1.1 in [10] are satisfied. Thus the scheme is strongly 
parabolic in the sense of Varah [9], and we shall rnaxe use of 
thi~ fact when treating the boundary conditions. 

We proceed by analyzing the scheme (2.13). We maxe the 
additional assumption that the eigenvalues of 



+ Bo >^«2 


+ Cw 


2 

2 



- 13 - 


are real. This assumption is valid for the Navler Stoxes equations, 
for example. Condition (ii) of theorem 1.1 in [10] is satisfied 
because nc w h is given by (3.8) with 


x = 2yAt (-^4 
(Ax) 2 
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and therefore, only one root lies on the unit circle. It remains 
to prove that the eigenvalues of the Fourier transform of (2.15) 
do not lie outside the unit circle. These eigenvalues are given, 
again, by (3.1) "ith 


-M 


r£r (i u + 




At 

Ax Ay 


8 < fi 2p>x< S 2 p y 


and 


X = 2vAt ( + 

(Ax ) 2 


py) 


9 


Note that ^p)x and <%>, are negative and (^ p ) x 

(£i ). are purely imaginary. For usual difference approxirr tic r.r 

y is positive. 

Condition (3»2) implies that v must be chosen such that 


( 3 - 10 ) 


At 

( a *y 


( D 2p^X 


At 


( Ay 


Ap 

(D 2 


2p^y ' 


At 

Ax Ay 


A] 

B (D^ 


At 

(D 


2p ; x v 2p y 


< 4yAt ( -8^4 + 

(Ax) 



) • 


By the parabolicity condition (2.1M) it is clear that (3.10) is 
satisfied provided that: 
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( 3 . 11 ) 


y > v. 



For the approximations (2.7) one can provq that (3.10) is satisfied for 


'■ 


77 

\%\<* 


■ 

2p v 


FJI 


A similar analysis holds for (2.18) and therefore we have completed 
the "tability proof of the generalized Du Fort-Frankel methods. 

We will new discuss briefly the effect of lower oi der terms 
in the eouation. The advantage of the Du Fort-Franxel scheme is 
that i* can be combined in a natural way with the Lea-;- Frog scheme 
if the equations jnder cons iderat ion have lower order terms. 

As an illustration we discuss the equation: 


(3.12) 


vl^ - Au 


+ <7U 


XX 


The scheme will be 


( 3 . 13 ) 


, n+1 n-1 

u ,i - u - 

U Urn , 

2 At 


, _1 n 

"V.l 


.2 n , .. 
" ?P U 3 ‘ Y<n '\l 


n+1 _ n n-1 , 

-2u 1 + u, ) 


The effect of the term AD^p on the amplification mat) ix will 
be 6i$f(5) in the upper left corner, where A is very small. 
Therefore, due to the discussion in [6, Sec. r.3] , the scheme 
’'•emai.os stable. However, if A depends on u and o is very 
small, we need a te^m that will maxe the scheme dissipative. 
Kreir.s and Oliger [^j suggested the form 


ORIGINAL PAGE IS 
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“ n+1 - (I + ^r(D + DJ 3 )u n_1 + 2 k aj3 c (i-K d J u '' • 
where 


n u * + i " ul J i 
D U*} = — itl ill 

o J 2 


A trivial calculation r>hows that a sufficient condition for 
stability is 



1 n~ 

2 V 0.Tf7ob‘ 


e < 1 . 


4 . The Fourier method for periodic boundary conditions . 

Let N be a natural number. Ax = - and = * 

j = 0,1, ...,2N . Consider a function u(x) such that u(x4-i) = u(x). 
An accurate method of approximating u at x = x. is to 

XX J 

interpolate u(x.) by the trigonometric polynomial 

•1 

N 

p A 

(4.1) u(x.) = \- «fi)cxp(:-ri«x .) , 


where 


u(®) - 


r~ 


A x . u ( x . )exp (-2iri®x . ) 

V 0 


and to differentiate this polynomial to get 


(4.2) 


* 2 u 


dx 


= V -4-rr^V '&(®)exp(2irimx .) 

J 


X=X .1 ®=-N 


This approximation can be achieved by two Fast Fourier Transforms 
and N complex multiplications. 



With 


- (u(x Q ufXgjj)) , 

a 2 «(x 0 ) i 2 n(^) 

v h - ( m \ji >•••> — rr— > 


bx 


C^X 


the above method can be written as = T^T\? N « S^u^. , where 
(4.3) A = diag(-4v 2 *T, -4* 2 (N-1) 2 , ... , -4 tt 2 ,0,-4v 2 , ...» -4ir 2 N 2 ) 


and 


(4.4) 


T m * ^Ax exp (-2vi(k-N)x t ) . 


It is obvious that T is a unitary matrix and therefore we have 


Lenaa 4.1 : is Hermitian and Its eigenvalues are all negative 

except one which is zero. 

Consider now the equation 

(4.5) “ t - <™*x > 

a 

u(x,0) = f(x) * T ^(:u)exp(2Tri®x) , 

®=-N 

u(0,t) = u(l,t) . 


V7e approximate (4.5) by 


(*. 6 ) 


-tf* 1 - -a"- 1 

Tt N _ e r * n M 2,~»n+1 „ n -»n-l . 

* 0S N U N * oyr (^k ' ?U N + U N > * 


2 At 

which can be written 
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(". 7 ) 


f *r 






where 


1 - V "! 

*Y T I 

I 

I 

C I 


The eigenvalue. • •: G N ^re given by equ&t-’ n (2.1) with 


2 2 

y = y(«) * 4Atotu n , x *= 2y<rAti; 


if now v > , then by Lemma 3.1 all the eigenvalues of 

except a simple one, lies inside the unit circle. In order to 
p-ove that 



( 4 . 9 ) 


II &S ii ^ K 


where K is not s function of either n c.«~ N , we shall prove 
that G can be diagonalized by a similarity transformation which 
is independent of N . Define first 

* TO 

(4.10) T = 

0 T 


Then 


ORIGINAL PAGE IS 
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A A 


TOT* - 



1— X 

1+x 


I 


A 

■ G . 


0 


The eigenvectors of 


A 


G 


are of the form 



0 

1 

0 


0 



0 


± 

X 2N+1 

0 


1 


where 



are solutions of (3.1) with 


y 


y(wi, * -N, -N+1,...,N. 


If 



t 


where 


* diag ( ^ 2 ***** i 1 ^ 


we have 
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~ (L + - L" ) _1 -L"(L + - L*) -1 

R" 1 « 

_-(L + - L")” 1 lV-L')’ 1 . 

|| R || and || R* 1 || are bounded Independently of N and R _1 GR 
Is diagonal* which completes the proof. 

The generalisation of this method to variable coefficients 
is trivial. Moreover, because of the parabolic ity there are no 
stability problems. This is not the case for hyperbolic equations, 
see [3]. 
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5. Boundary Conditions 

In this section we will treat only 4th order approximations 
to the scalar equation (2.1). The results can also be applied 
to systems, but the boundary conditions must then be stated in 
terms of those variables corresponding to the diagonalized 
system (If such a form exists). We will always consider the 
quarter space problem 0 < x < •, t > 0 and the theory by 
Varah [9] will be used. We begin with the Dlrlchlet boundary 
condition 

u(0,t) - g(t) 

The stability proof will be applied to the more general 
scheme 

(5.1) l (a. I ♦ B,.D u 2 )u, n “ k - 0 

k«-l K J 

p 

where is one of the 4th order accurate operators treated 

in Sec. 2. (5.1) is assumed to be consistent with (2.1) and 

staole for the initial value problem. Furthermore we assume 
strong parabolicity , which was shown to be fulfilled by the 
Du Fort-Frankel scheme in Sec. 3. For the explicit operator 
defined by (2.7b), the numerical boundary conditions 
will be 

(5.2) u Q n - g n 

k L [ V*W- (I - r? D . 2 >V k - 0 


(5.3) 
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In (5.3) the centered difference operator used for inner 
points Is substituted by a non-centered 3rd order accurate 
operator. This lower order accuracy at the point x » Ax 
should not effect the overall accuracy, (see Sec. 6). 
Connected with (5.1) is the resolvent equation 

(5. ft) T (z )5j - D + D_ ( I - D + D_)Uj - 0, j - 2,3,... 

where 

*c*> - ! v r -/f v r ' k - i=i > i. 

k--l K / k*-l 

and the characteristic equation 
(5.5) t(z) - f(x) - 0 


where 

f(o - ( - K ~^ - 2 (1 


(K-iy 

"T5F" 


) 


(5*5) has for |z| >1 two roots ic^, < 2 inside the unit 
circle (see [9]). We can now prove 


Theorem 5.1 . Assume that at least one of the roots, < 2 (z) 
say, satisfies ** IC 2^ Z ^ ** *0 for l z ! - where 

Kq is the root to 

(5.6) ft/ - 9< 3 - 5* 2 - 15 k ♦ 1 - 0, 

which is Inside the unit circle. Then (5*1) with D^ 2 = 
is stable with boundary conditions (5.2), (5*3). 
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Proof * 

The homogenous boundary conditions for Qj are 

(5.7) u Q - 0 

(5.8) t(z)u 1 ♦ D + D_(I - 17 V> 5 1 * 0 

Looking for non trivial solutions In t 2 (0,«) (i.e., 

£ |u. | 2 Ax < ®) for |z | > 1, we write u. in the form 

J»0 J J 

(5.9) Uj « a(K^ - «c 2 ^), < 1 f* < 2 

where are defined by (5.5). 


The condition for the existence of non trivial solutions is 

obtained by inserting (5-9) into (5*8): 

9 U,-l ) 2 p (K -i ) 2 

(5.10) t(z)(k 1 -ic 2 ) - [(iCj-D^d y 2 ) " (< 2 ~ 1 ^ ^ H ^ 


= r> 


By combining (5.10) with (5*5) it can be shown that ic^ 
must fulfill 

(k-,- 1)' (iCp-l)-* 

(5.11) — l ^ * 0 

*1 *2 

together with one of the equations 

(k -l) 2 (Kp-1) 2 

(5.12) — i + — £ 12 


or 


(k^-1) 2 (k 2 -1) 2 


0 


(5.13) 
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We begin with the system (5.11), (5.12). By defining 

*1“^ c 

c * — - _ y we get < 1 * c k 2 > and from (5-11) 

* -c(c+l)(c 2 +l). 

(5.12) then gives the final polynomial for c: 

c 10 c 9 + 2c 8 + 3c 7 ♦ 16c 6 + 3c 5 + 16c 11 + 3c 3 + 2c 2 + c ♦ 1 » 0 

The roots were obtained by a computer program. They are all 
such that at least one corresponding is outside the unit 

circle, which is a contradiction 


An easy calculation shows that the system (5.11), (5.13) 
implies that and in this case the form ( 5 . 9 ) of 

the solution does not hold. For double r^.ots of (5-5) the 

a -I 

form is u^ * aj< . By inserting it into (5.8) and using 
(5.5), we obtain immediately the condition 

4e 5 - 13*^ + 4 k 3 - 10ic 2 + 16 k - 1 = 0 

for a nontrivial solution. k=1 is one root, but is ruled cut by our 
assumption. The remaining deflated polynomial is (5.6), and the 
theorem is proved. (It should be noted that the assumptions m r • 
theorem could be weakened, since Varah's stability condit.^n perm t s 
non-trivial solutions of a certain type.) 


ORIGINAL PAGE IS 
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Corollary . The Du Fort-Frankel scheme is stable with boundary 
conditions (5-2), (5*3). 

Proof . ic»l is actually a double root of (5.5) for z=l but 
a perturbation calculation shows that only one of them is 
inside the unit circle for |z| > 1. Therefore the first 
assumption of the theorem is fulfilled. 

The root »c Q of (5.6) with |k q | < 1 is real and 
located in the interval [0.06, 0.07]. This corresponds 
to a value of f(ic) which is less than the stability limit 
16/3, and, therefore |z| <1 (see Sec. 3.). 

It is easily verified that the assumptions in the theorem 
are fulfilled by, for example, the following schemes 

u n+1 - wQl u n , 

y Ql u n+1 - u n , 

(I - euQ-^u 11 * 1 « (I + (l-0)yQ 1 )u n , 

((1+6)1 - vQ 1 )u n+1 + (1-20 )u n + 0U 11 ' 1 = 0, 

where it is assumed that u « JAt - a - and 0 are In the 

(Ax) 2 

stability intervals for the initial value problem. The 
theorem can also be generalized to several space dimensions 
in the sense that after a Fourier transformation over all 
space variables except x, the stability conditions of Varah 
[9] still are fulfilled uniformly in the dual variables 



25 


C 2 » . This depends on the fact that the drop 

out of our calculations as z and u do, and the final 
polynomial will be independent of all parameters. 


Next we will treat the operator Q 2 defined by (2.8) 

One boundary condition for the difference scheme obviously 
must be (5.2). This condition will also be sufficient to 
define the solution if we multiply (5.1) by D + D _) 

and solve for u n+1 directly. In this case stability 
follows immediately, since the solution to the resolvent 
equation has the form Uj « aic J , and it cannot satisfy u Q « 0 
if a ^ 0. However, this procedure will become inconvenient 
for a real problem, where the coefficients depend on x, t 
and very likely even on u. 


In a practical application for an explicit scheme, the 
second derivative (Ax)~^Sj n is computed from u^ n by 
solving the tridiagonal system, 

(I + D + DjS. n - D + D_ Uj n , J-1,2,... 

For an implicit scheme, the system 


(I ♦ D + DJ Sj 


n+1 - D + Djlj n+1 » 0 


n+1 . ft o n+1 
a 0 u j 6 0 S J 


-l (<Vj n ‘ k + 


k«0 


S n " k } 
k J } 


J-1,2, 
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is solved for Uj n+1 , Sj n+1 . In both cases, a boundary 
condition is needed for S Q . One way of proviiing that is 
to express S Q in terms of Uj , and this is also the simplest 
way for an explicit scheme. 

By using a 3rd order accurate one sided formula, w<_ get 
the boundary condition 

(5.14) S Q n - (35u 0 n - 104 U;L n + ll4u 2 n - 56u 3 " + llu M n )/12. 

2 

Theorem 5.2 . The scheme (5*1) with = Q 2 is stable with 

the boundary conditions (5.2), (5.14). 


Proof . 

The resolvent equations for Uj , Sj are 


T(Z)S J - Sj 

(I + D + D_)Sj * D + D_Uj j 


u*l,2. 


with the homogenous boundary conditions 

S o*°. S o - J 0 e J G j* 

e^ being defined by (5.14). After eliminating Sj , these 
equations can also be written 

(5.15) t(I ♦ D^DJuj * D + D -“j’ J-2,3,... 



( 5 . 16 ) 
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(5.1?) tu x ♦ j^{tu 

The solution in 
the form 


2 - 2t “i ♦ ' “a - 2i 


(0,«) to this system has for 


1 


1*1 > 1 


(5.18) Uj • a(*^ - 6^) (fi^ Is the Kronecker symbol) 
where k satisfies 

2 

(5.19) t(s)(k ♦ ^fp-) - (*-l) 2 - 0 

When using (5.18) , (5.19) in (5.17) we obtain the condition 
for a non trivial solution: 

(k 2 ♦ 10w ♦ 1) T e-K^** 1 ♦ 14X» - 0 
J-l 3 

For our choice of this equation Is 

(5. 20) 11k 5 ♦ - 435k 3 ♦ 980k 2 - 926k ♦ 162*1 » 0 


which has no root inside the unit circle, and the theorem 
is proved. 

We consider next the boundary condition 
u x (0,t) ♦ bu(O.t) » g(t). 

For both of the operators and Q u x (0) * s 

approximated by a one-sided 4th order accurate formula 

(5.21) (-25u q ♦ 48u x - 36u ? ♦ . >u 3 - 3uu)/12Ax ♦ bu Q « g. 

In connection with the second boundary condition wi^.1 
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be (5*3)* The stability of the scheme has not been verified 
theoretically. However, numerical experiments show no 
evidence of instability (see Sec. 6.). 

For the implicit operator Q^, the second boundary 
condition will be (5.14), and in this case we can prove 
the following 

Theorem 5. 3 . 

2 

The scheme (5.1) with = Q 2 is stable with the 

boundary conditions (5.1*0, (5.21). 

Proof . 

The proof follows the same lines as the proof of 
Theorem 5.2, and we do not give the details here. Keeping 
the parameter a defined in (5.13), the final equation 
corresponding to (5*20) now is 

(5.22) P(ic)a = (602k 5 + 2376k 1 * - 24064k 3 ♦ 56764k 2 - 71546k ♦ 35363)a 

» 0 , 

which has 4 roots outside the unit circle, and one root 
at <*1. By the assumption of strong parabolicity , the 
corresponding value of z obtained from (5.13) is 1, and 
for this z-value, k» 1 is a double root of (5.19). 

Therefore k ■ 1 ♦ d ( / | z-1 j ) . The stability condition by 
Varah is that the parameter a can be estimated from 
P(Oa * g by f a | < | g | / /| z-1 1 . This estimate is valid 
for our case since k*1 is a simple root of P(»c). 
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6. Numerical Experiments 

The generalised Du Port-Prankel scheme was tested for 
the Burgers equation 

(6.1) u^ ♦ (u-l/2)u x - Jtt xx , “5 < x < 5, 0 < t, 

using the two 4th order accurate operators Q^, Q 2 defined 
In Sec. 2. The Initial function was 

u(x,0) - 1 - (x-5)/10, 

and the boundary conditions were 

(6.2) u(-5,t) - 1 

(6.3) u(5,t) - 0 

The problem (6.1), (6.2), (6.3) has the steady state solution 

(6.4) v(x ) » j<1 - tarh — - ). 

The time Integration was stopped when the condition 

(6.3) max |u. n+1 - u, n j < At • 10“ 5 
i J J 

was fulfilled. The error c * max|uj n - v(Xj)| In listed 

J 

in table 6.1. Parameter values used were o * 1/8, 

Ax ■ At *» 0.2, y ■ 4/3. € denotes the time when (6.5) 

was first satisfied. 
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t 10.6 lO" 4 3.4 1<T 4 


t 26.6 26.4 

A 

Table 6.1 The error e and steady state time t for the 
different operators. 

Our scheme was also run with ax * At * 0.1, and the 
error was found to be 16 times smaller in accordance with 
the 4th order accuracy. 

With the same initial function but with the boundary 
condition u x (-5) * 0 instead of (6.2), we obtained the 
steady state solution u(x,®) = 0 for both operators Q^, 
Q^. In both cases the scheme was not showing any signs of 
instability, which possibly could have occurred from the 
boundary condition at x * -ii, when usinp - the operator 
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